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ABSTRACT 

We consider the response of a protostellar disc to a tidally induced warp and the 
resultant changes in the spectral energy distribution (SED). We argue that for typical 
protostellar disc parameters the warp is communicated through the disc in a wave-like 
fashion. We find that the main effects of the warp tend to be at large radii (R > 30 
AU) and, for sufficiently small viscosity, can be quite long-lived. This can result in 
non-uniform illumination of the disc at these radii and can induce significant changes 
to the SED at wavelengths A > 100/im. 
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1 INTRODUCTION 

Protostars have circumstellar discs from which they accrete 
matter. Once the local accretion rate onto the protostar falls 
below a certain value, given approximately by 

where M is the stellar mass, R, the stellar radius and L* 
the stellar luminosity, then most of the energy radiated from 
the disc comes from re-radiated stellar emission, rather than 
its own internally generated accretion energy. Such discs are 
known as passive discs and most observed protostellar discs 
fall into this category. The nature of the spectrum emitted 
by the disc depends on how much of the stellar flux can be 
intercepted and at what radii this occurs. A flat, infinitesi- 
mally thin disc, has at large radii a temperature distribution 
of the form T<j oc i? _3//4 which gives rise to an SED of the 
form \F\ tx A~ 4 / 3 . Almost all protostellar SEDs are flatter 
than this, implying that real discs intercept substantially 
more of the stellar flux. For this reason, Kenyon & Hartmann 
(1987; see also the review by Dullemond et al., 2007) argued 
that protostellar discs are substantially flared, that is, that 
the disc scaleheights H are such that d(H/R)/dR > 0. Since 
then more sophisticated models of the flaring have been in- 
troduced (Chiang & Goldreich, 1997; D Alessio et al., 1999; 
Chiang et al., 2001; Dullemond & Dominik, 2004a,b), to- 
gether with truncated discs which display a heated inner 
rim (Natta et al, 2001; Dullemond et al. 2001; Muzerolle, et 
al., 2003; Isella & Natta, 2005). These models are generally 
successful in providing explanations of the observed SEDs 
in the near and mid-infrared. 



These models, however, all assume the disc to be cylin- 
drically symmetric, and this seems to be a reasonable as- 
sumption in the inner regions (Monnier et al., 2006). Stars 
typically form within clusters, and so even single stars are 
likely to have undergone close encounters with neighbour- 
ing stars at a young age. Very close encounters are likely 
to be dissipative and lead to disc truncation and perhaps 
even capture and binary star formation (Clarke & Pringle 
1991, 1993; Hall, Clarke & Pringle, 1996; Boffin et al., 1998; 
Watkins et al., 1998a,b; Pfalzner et al., 2005; Moeckel & 
Bally, 2006). More distant encounters occur more frequently, 
but can still provide the disc with significant warping. We 
note that there is observational evidence for such disc distur- 
bances (Quillen, 2006; Lin et al., 2006; Cabrit et al., 2006) 
and for the warping of the disc at large radius (Hughes et 
al., 2009). 

Viscosity in the outer regions of protostellar discs is 
likely to be small, and typical figures of the dimensionless 
viscosity parameter a (Shakura & Sunyaev, 1973) in the 
range 10~ 2 — 10 -4 are thought likely, e.g. Terquem, 2008. 
This is less than the disc thickness ratio H/R, and in this 
case such warps are propagated as waves (Papaloizou & Lin, 
1995). With such low viscosities, these waves can exist for a 
significant fraction of the protostellar lifetime (Section [2]). 

In this paper, we investigate the effects of these waves 
in a very simplified manner. We do not attempt to model 
perturbations to the structure or SEDs of discs with flares 
or other configurations. Rather we consider the effect on the 
shape and on the resultant SEDs emitted by a disc which 
is initially thin and flat. Although this is not likely to en- 
able us to provide realistic fits to any observations, it does 
enable us to see at what wavelengths any complications de- 
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riving from such perturbations might manifest themselves, 
and to provide an estimate of their relative magnitude. We 
note that the effect of warping on the SED of a protostellar 
disc has been considered previously for the case of a disc in 
a misaligned binary system by Terquem & Bertout (1993, 
1996), but that they do not consider the dynamic response 
of the disc, nor its subsequent evolution towards alignment 
(e.g. Bate et al, 2000). 

In Section [21 we provide estimates of the periods and 
lifetimes of warp-like perturbations in a protostellar disc. 
In Section [3] we compute the time-dependent behaviour of 
warps which can be generated using a simple model of the 
flyby of a perturbing star. In Section [4] we compute the 
temperature distributions and resulting SEDs for these disc 
models. We discuss the implications of our results in Sec- 
tion [5] 



2 WARP PROPAGATION AND LIFETIME 
ESTIMATES 

We consider a protostellar disc in orbit around a star of mass 
M. For a Keplerian disc the angular velocity of disc material 
is D. = (GM/R 3 ) 1/2 . We take the disc semi-thickness to be 
H, so that at radius R 



(2) 



H = c ± _ 
R ~ Rfl' 

where c a (R) is the local disc sound speed (Pringle, 1981). 
Typical values of H/R are expected to be around 0.05 - 
0.1 (Bell et al., 1997; Terquem, 2008), independent of radius 
R. Typical estimates for the dimensionless viscosity param- 
eter a for such discs are in the range a ~ 10~ 2 — 10 -4 
(Hartmann 2008). Thus for protostellar discs we expect that 
a < H/R, and therefore that warps propagate through them 
in a wave-like manner (Papaloizou & Lin, 1995). The prop- 
agation speed for these waves is (Papaloizou & Lin, 1995; 
Pringle, 1999; Lubow, Ogilvie & Pringle, 2002) 

Vw = ^c s . (3) 

Suppose that a warp is induced in the outer disc regions 
by the fly-by of some other stellar object. If the relative ve- 
locity of the encounter is V and the distance of closest ap- 
proach is D then the characteristic frequency of the induced 
warp-like disturbance is 

V 

"» ~ (4) 
and the wavelength of the disturbance is 



A w (i?) ~ — 



c s (R) 
2V 



D. 



(5) 



We expect the warp to propagate inwards to a radius 
7? cr it at which R — X(R). At that radius the propagating 
warp wave is reflected, and within that radius the disc is 
flat, but tilted by the warp. This is because the wavelength 
of the warp is greater than the size of the disc in the inner 
regions and therefore cannot communicate the warp. How- 
ever angular momentum is still communicated into the inner 
disc, so the disc tilts without warping. If H/R « const., in- 
dependent of radius, then we find that at radius -Rcrit 



f^(-Rcrit) ~ "77"^w 
tl 



(6) 



(7) 



If we make the further assumption that the fly-by is 
parabolic, so that to a first approximation 

then we find that 

Thus we expect that while the inner disc regions R < _R C rit 
will have a time- varying tilt, it is only the outer disc regions 
R > -Rcrit that will show a significant warp. As a numerical 
example, if we assume that the fly-by passes the disc at a 
distance of D — 250 AU, and H/R — 0.1 then we have 
.Rcrit « 34 AU. 

The time taken for the warp to propagate to the disc 
centre and back is 

4J? ou t 



(9) 



where Rout is the radius of the outer disc edge. Equivalently 
this may be written as 

2 7? 



tt'H 



■ ■Pqu 



(10) 



where P out is the orbital time at the outer disc edge. Thus 
for an outer disc edge at 100 A.U., and 

H/ Re* 0.1, (11) 

the warp crossing timescale is around 10 4 yr. 

The decay timescale for the warp is given approximately 
by (Lubow et al., 2002) 



fda 



Pout 

2ira 



(12) 



Thus for a disc with R out = 100 AU, we would expect the 
warp to last for a time of around 1.6 x 10 6 (a/10~ 4 )" 1 yr. 



3 WARP GENERATION - NUMERICAL 
RESULTS 

3.1 The warp equations 

We model the generation and propagation of the warp using 
the linearised equations derived by Lubow & Ogilvie (2000) . 
These describe the warp in terms of the local unit tilt vector 
l(R, t) which is a function of radius and time. For a Keplerian 
disc these are 

19G + T, (13) 



dt R dR 



and 
<9G 



+ qJ]G 
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(14) 



dt 4 dR' 

Here G is the internal torque which acts to realign tilted disc 
annuli. T is the external torque acting on the disc which in 
this case is caused by the stellar fly-by. In addition S is the 
disc surface density, and P = J pdz = Q. 2 T,H 2 is the verti- 
cally integrated pressure. We note that these equations do 
not make allowance for dispersive or non-linear effects in the 
wave propagation (Ogilvie 2006). For the waves we discuss 
here dispersive effects may be of marginal significance, but 
only at the outer disc edge. However for the large ampli- 
tudes required to substantially change the SED, non-linear 



Warps in protostellar discs 3 



effects may well play a role. However, to take these effects 
properly into account would require a full hydrodynamical 
model of the disc, which is beyond the scope of the current 
investigation. 

To solve these equations numerically we adopt the ap- 
proach of Lubow et al. (2002). The disc is taken to lie in the 
XY-plane, so that for small tilts we have to first order 



E(R) oc R~ 



(22) 



1=(J.,I»,1), 



(15) 



ox, i-y ^ j-- The torques G and T lie in the XY— 
plane. We can then use W = l x + il v , G = G x + iG y and 
T = T x + iT y as complex representations of 1, G and T, 
respectively. We then replace the quantities W and G by 
the quantities A(R, i) and P(R,t) defined by 



W = 
and 



D* 



G = -YjH 2 R 2 Q.A* , 

where the asterisk denotes a complex conjugate. 
Equations [13] and [14] then become 



(16) 



(17) 



dP 
~dt 
and 
dA 
~dt 



Si 
2 



d 



Ei? 1 / 2 // 2 OR 



(T.R 



1/2 



H 2 A) 



-aQ,A 



(19) 



1 dP _ P 
2~dR ~~ ~R' 

We truncate the disc at R = Rout and do not allow for 
dissipation at the outer boundary. This is because the wave- 
length of the warp is long in comparison to any smoothing 
of the outer edge and so the warp will only see it as a sharp 
edge and reflect back inwards. 

3.2 External torque due to fly-by 

We now obtain an expression for the torque density T{R,i) 
exerted upon an annulus of the disc during a stellar fiyby. 
The torque density exerted on the annulus at radius R with 
tilt vector 1 by a mass of Mi at position vector Rb is given 
to first order in R/R h by (Lubow & Ogilvie, 2000) 



GM 2 
2Rt 



3 ' 2 \R h J 



E_R(R b -l)(R bA l), 



(20) 



where b\L is the Laplace coefficient. We approximate this 
in the form 



{Rl + R?f/- 



■(R b -l)(R b Al). 



(21) 



We note that since we are using the approximation of only 
considering linear warp waves, it is only the T x and T y com- 
ponents which are relevant for the computations here. 

3.3 The disc shape and its evolution 

We start with a disc which has an outer radius of Rout = 100 
AU. We take the inner radius to be at _Ri n = 0.01AU, and 
note that this is much less than the expected value of -R cr it 
(equation lll|l within which we expect the disc to be tilted 
but not warped. 

We take the surface density of the disc to be a power-law 
of the form 



(c.f. Andrews et al., 2009; fsella, Carpenter & Sargent, 2009) 
and take H/R = 0.1 = const. This then implies that the 
sound speed is a power-law of the form: 

-1/2 



oc R~ 



(23) 



We apply a small amount of damping by taking a — 10~ . 
We note that for these disc properties, we expect the warp 
amplitude to vary in such a manner as to keep the flux (-Fx) 
of angular momentum associated with the warp constant, 
where 



F± oc \W\(GMR) 1/2 T,R 



(24) 



For the disc properties we have chosen this corresponds to 
\W\ « const., so that the vertical local disc displacement 
due to the warp is roughly proportional to radius. 

The disc evolution equations [18] and [19] are evolved nu- 
merically using a leapfrog scheme. Boundary conditions are 
chosen to ensure that the inner and outer edges of the disc 
are stress-free. 

To illustrate the possibilities we consider two idealised 
flybys. We take the orbit of the perturber to lie along a 
straight line, constant velocity path. In practice, of course, 
the interaction is more dynamical than this (e.g. Cabrit et 
al., 2006; Moekel & Bally, 2006), however we can, to first 
order, obtain the correct warping structure for the disc. Thus 
in both cases we take the path to be of the form 



R b =a + V(t-t ). 



(25) 



Here to is some arbitrary time at which the perturber is at 
position vector a. 

3.3.1 Model 1 

For Model 1 we take the fiyby to lie along a straight line 
perpendicular to the disc of the form 



R b = (300, 0, V(t - i ))AU. 



(26) 



Here to is the time of closest approach when the perturber 
is a distance of a = |a| = 300 AU from the disc centre, and 
200 AU from the disc edge. As can be seen the velocity V is 
in the Z-direction. We take its magnitude V to correspond 
to the escape velocity from the point of nearest approach, 
so that 

1/2 



T , f2GM\ 1 / 2 ( M\ ' ( a \-V 2 

V=( =2.44 —— — kms \ 27) 

V a J \M a V300AU/ v ; 



Because the path of the perturber lies solely in the XZ- 
plane, the disc annuli are tilted only about the F-axis. Thus 
we can envisage the disc warp by considering a cut through 
the disc in the X Z-p\&nc. 

The results for Model la for which we take the mass 
of the perturber to be M2 = IMq are shown in Figure [T] 
At time t = the disc is assumed unperturbed and the 
perturber has a z-coordinate of z = —1000 AU, and thus 
to ~ 2000 yr. The outer edge of the disc is initially pulled 
downwards. Then as the perturber passes the disc (passing 
through the initial disc plane, at time t = to ~ 2000 yr) the 
outer edge of the disc starts to be pulled upwards. Thus the 
disc is set oscillating with an amplitude at the outer edge of 
around 5 AU. The period of the oscillation is 
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P = 2 



dR 



4 R 



c s (R) 3c s (R a 



(28) 



t-1500yrs 



For the disc considered here this gives P 2100yrs. 

For Model lb we explore a more substantial warping of 
the disc. For this we use a = |a| = 200 AU and we take 
the mass of the perturber M2 = IOMq. We note that this 
changes the velocity of the flyby to V = 9.5 km s _1 . As can 
be seen in Figure [2] this produces an amplitude for the warp 
wave of approximately 35 AU. 

3.3.2 Model 2 

In Model 2 we take the path to lie in the direction V = 
(1, 1, 1), that is, at an angle to the disc plane. The perturber 
is initially at position 



R b (t = 0) = (-700, -1000, -1000)AU. 



(29) 



The perturber is then given a velocity of V = |V| = 2.7 km 
s _1 which corresponds to the escape velocity at the point 
of closest approach R b = (200, -100, -100) AU. Note that 
this path passes through the XY-plane at the same point as 
the path in Model la. Because this path is not in the XZ- 
plane, the torque acting on disc annuli is now no longer 
simply in the Y-direction. Thus the disc is not only tilted 
but also acquires a twist. For comparison with Model la we 
take M 2 = 1M . 

As the disc is twisted we can no longer view it via a 
cut through the XZ-plane. To illustrate this, in Figure [3] 
we show the disc in three dimensions, for Model lb, and 
compare it to Figure [U which is a typical shape of the disc 
for Model 2. The disc shapes for Model la and lb are the 
same but with differing amplitudes, so we use Model lb here 
to illustrate the disc shape more easily. 



4 DISC HEATING AND RESULTANT 

SPECTRAL ENERGY DISTRIBUTIONS 

4.1 Temperature profile 

In order to keep things simple we assume that the disc is 
thin, and therefore ignore any flaring. In this way we ensure 
that any re-emission of radiation from the disc at large radii 
(long wavelengths) is caused predominantly by the warp. 
As we have noted, the inner disc (well within the radius 
R cr it ~ 30 AU) remains flat, but with variable tilt. Thus for 
the inner disc we may take its surface effective temperature 
to be T d (R) where 



Ti(R) = Tt 



Hi) 



~~R 



ffi 



(30) 



where T* is the stellar effective temperature and R, the 
stellar radius (Friedjung, 1985; Kenyon & Hartmann, 1987). 

For the warped part of the disc, we write the local unit 
vector representing the disc tilt as 



l(i?) = (cos 7 sin /?, sin 7 sin j3, cos (3) , 



(31) 



so that the angle /3(R) represents the tilt of the disc normal, 
relative to the Z-axis, and the angle 'y(R) represents the 
azimuth of the tilt. Since the warp starts only at a distance 
R > -Rent 3> R* we are able to treat the illumination of the 
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Figure 1. A sample of the XZ-profiles of the disc at various 
times for Model la, showing the induced warp and its propagation 
through the disc. For time f = 0, the disc is entirely flat in the 
XY— plane. Here we only show the positive X-axis. In the XZ- 
plane, the disc is anti-symmetric around x = 0. The axes are in 
AU. Note that the inner region of the disc (R < 30 AU) is tilted 
but not warped. 
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Figure 2. An example of the XZ-profile of the disc for Model 
lb, showing that in this case the warp is more substantial at the 
outer edge of the disc. The axes are in AU. 
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Figure 5. Temperature profile for the disc in Figure [3] from 
Model lb. The axes are in AU, with the temperature coloured by 
the scale shown in Kelvin. The central star (with T* = 5000K) 
and innermost parts of the disc are omitted to give a higher res- 
olution over the disc. Note that only one half of the disc appears 
illuminated to an external observer. 




Figure 3. We show here a representative disc shape in three 
dimensions for Model lb. The image was generated using the 
graphics package POVRAY. The top image is the XZ profile of 
the 100 AU disc, viewed and illuminated from (0, -200, 0). The 
second image is the same disc but viewed from (0, —200, 50), and 
illuminated from (0,5,0). Because the path of the perturber lies 
in the XZ— plane, each disc annulus is tilted about the Y— axis. 
Thus the disc is warped but not twisted. 




warped part of the disc by the central star as coming from a 
point source. In this case (cf. Pringle, 1996) the disc surface 
temperature can be written as 



n 



| sin <f>R — cos <f> sin ft R^ \ 



1+ (cos <t> sin P R§± - sin 0-R§§) 2 



1/2' 



(32) 



This is valid only if self-shadowing of the disc can be ig- 
nored. In practice this needs to be taken into account. Thus 
this temperature is applied only to those parts of the disc 
which are tilted towards the central star and which have 
a view of the central star which is unobscured by other 
(warped) parts of the disc. Parts of the disc which are ob- 
scured are assumed to have negligible temperature. This 
computation is carried out in a straightforward manner as 
in Pringle (1997). 

Once the temperature distribution of the disc is found, 
we compute the emitted disc spectrum simply by assuming 
that each disc element emits locally as a black body. For 
the SED we assume an observer placed on the Z-axis for 
simplicity. If the disc were to be viewed from an angle there 
would be an apparent reduction in the magnitude of the 
flux emitted by the disc. If that angle is increased enough, 
then the central star and large portions of the disc would be 
obscured from the observer by the warp. 




Figure 4. We show here a representative disc shape in three di- 
mensions for Model 2. The image was generated using the graph- 
ics package POVRAY. The top image is the XZ profile of the 100 
AU disc, viewed and illuminated from (0,-200,0). The second 
image is the same disc but viewed from (0,-200,50), and illu- 
minated from (0,5,0). Because the path of the perturber is not 
perpendicular to the unperturbed disc, the individual disc annuli 
are tilted about different axes in the XY— plane. This means that 
the shape of the disc is twisted about the Z-axis. The central 
region is tilted in one plane and the outer rings are slowly twisted 
away from this plane, producing the shadowing effects. 



4.1.1 Model 1 

In Figure [5] we show a representative temperature distribu- 
tion found for Model lb, and in Figure [6] the resultant spec- 
tral energy distribution (as viewed by a distant observer in 
the direction of the Z-axis) from the disc structure shown in 
Figure [3] For comparison we also plot, on figure [(J an SED 
from Model la and the contribution to the SED from the 
central star. Because of the induced warp, the SED displays 
a secondary peak at around 100 /J,m. For shorter wavelengths 
the radiation comes from the stellar surface and the flat illu- 
minated part of the inner disc. The magnitude of the warp is 
larger in Model lb, creating a more substantial flux from the 
warped disc at around 100 /im as opposed to the weaker flux 
from Model la. In Model 1 all disc annuli are tilted about 
the same axis, therefore outside the central region (taken to 
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Figure 6. SED in terms of \F\. The solid line is for the disc in 
Figure [3] from Model lb. The dot-dashed line is a typical SED 
from Model la, and the dashed line represents the star which is 
taken to be a black body with T» = 5000 K. The contribution 
from the star peaks at about 0.6 (im. The inner disc has a similar 
maximum temperature and a log wavelength dependence of the 
form \F\ oc A -4 / 3 . The contribution from the upturned part 
of the outer disc peaks at around 100 fim. As the warp is more 
substantial in Model lb, the flux in this peak is greater than the 
flux in the corresponding peak for Model la. As is to be expected, 
both disc peaks have less magnitude than the primary peak, with 
Model lb having a disc peak flux that is around 5 times less than 
the primary peak flux. 

be at the reflection radius), only one half of the disc appears 
illuminated. 

In general the disc illumination can be more compli- 
cated than this because of self-shadowing, and we illustrate 
some possibilities in Figure [7] 

4.1.2 Model 2 

As the flyby in Model 2 is misaligned with the original plane 
of the disc, the resultant disc is both tilted and twisted. 
This can be seen in Figure [8] which shows a typical temper- 
ature profile for Model 2. This exhibits a spiral pattern in 
the illumination profile. The resultant SED, Figure^ is not 
changed greatly by this, which is to be expected as the disc 
intercepts a similar amount of stellar flux when both tilted 
and twisted as opposed to just being tilted. The contribution 
of the heated outer disc again peaks at around 100 fim. 



5 DISCUSSION 

We have considered the response of a protostellar disc to a 
tidally induced warp, and the resultant change in the spec- 
tral energy distribution (SED). We have used a simplified 
analysis in order to emphasise the relevant implications. 

For the low viscosities currently envisaged for the outer 
disc regions we note that a tidally induced warp is propa- 
gated as a wave, and that such a wave can be comparatively 
long lasting (in our case more than a Myr for a = 10 -4 ). 
The warp is largest in the outer regions, and within some 
critical radius (in our case i? cr it ~ 30 AU) the disc remains 
flat but has a tilt which is variable on a timescale of order 
the propagation timescale of the warp (here a few thousand 




Figure 7. Some examples of the complicated illumination pro- 
files that arise due to the self-shadowing as the warp propagates 
through the disc. Note that the characteristic property of an il- 
luminated warp is that the disc surface brightness is lop-sided 
with respect to the central star. The images were generated by 
the graphics package POVRAY, with the source of illumination 
on the central stars surface. The discs here have a radius of 100 
AU. 



years). Note that any collimated outflow, or jet, originating 
from the inner disc regions might be expected therefore to 
vary direction on this timescale. 

In reality the passive disc that we are considering here 
will be flared. Flaring of the disc contributes to the SED at 
wavelengths longer than a few /jm. In this paper, in order 



Temperature Profile/K 




Figure 8. Temperature profile for the disc in Figure [4] from 
Model 2. The axes are in AU, with the colour scale giving the 
temperature in Kelvin. The central star (with T„ = 5000K) and 
innermost parts of the disc are omitted to give a higher resolution 
over the disc. Note that when the path of the perturber is not 
perpendicular to the unperturbed disc plane, the disc is twisted 
and so the pattern of illumination displays an apparent spiral 
pattern. 
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Figure 9. SED in terms of \F\ for the disc in Figure [4] from 
Model 2. 



to emphasise the effects of the warp, we have assumed the 
disc to be locally thin when computing the SED. In practice 
the two effects, that of the flare and that of the warp, must 
be combined. 

From our computations, in agreement with Terquem & 
Bertout (1993, 1996), we conclude that the warp is only 
likely to affect the disc SED at wavelengths of around 100 
/im and longer, but that the effects can be substantial. We 
also note that warp-induced variations in disc illumination 
have the effect of making the disc appear asymmetric with 
respect to the central star at radii R > i? C rit- Further, a twist 
in the disc, induced by a skewed flyby, can have the effect of 
introducing a spiral component into the disc illumination. 
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